Astronomical Journal, December 2002, in press 



Hubble Space Telescope Evidence for an Intermediate-Mass Black Hole in the 

Globular Cluster M15 — 
II. Kinematical Analysis and Dynamical Modeling^ 

Joris Gerssen, Roeland P. van der Marel 
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218 

Karl Gebhardt 

Astronomy Department, University of Texas at Austin, Mail Code C1400, Austin, TX 78712 

Puragra Guhathakurta, Ruth C. Peterson^ 
UCO/Lick Observatory, Department of Astronomy and Astrophysics, University of California at 
Santa Cruz, 1156 High Street, Santa Cruz, CA 95064 

Carlton Pryor 

Department of Physics and Astronomy, Rutgers University, 136 Frelinghuysen Road, Piscataway, 

NJ 08854-8019 



^Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope 
Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA 
contract NAS 5-26555. These observations are associated with proposal #8262. 

^Also at: Astrophysical Advances, Palo Alto, CA 94301. 



- 2 - 



ABSTRACT 

We have used the Space Telescope Imaging Spectrograph (STIS) on the Hubble 
Space Telescope (HST) to obtain high spatial resolution spectroscopy of the central 
region of the dense globular cluster M15. The observational strategy and data reduction 
were described in Paper I (van der Marel et al. 2002). Here we analyze the extracted 
spectra with a cross-correlation technique to determine the line-of-sight velocities of 
individual stars. Our final STIS velocity sample contains 64 stars, two-thirds of which 
have their velocity measured for the first time. The new data set triples the number 
of stars with measured velocities in the central projected R < 1" of M15 and doubles 
the number in the central R < 2" . We combine our data with existing ground-based 
data to obtain non-parametric estimates of the radial profiles of the projected rotation 
velocity, velocity dispersion, and RMS velocity ctrms- The results differ from earlier 
work in the central few arcsec in that we find that urms rises to ~ 14 km s"-*^, somewhat 
higher than the values of 10-12 kms"^ inferred previously from ground-based data. 

To interpret the results we construct dynamical models based on the Jeans equation 
for a spherical system. If the velocity distribution is isotropic, then M15 must have a 
central concentration of non-luminous material. If this is due to a single black hole, 
then a fit to the full velocity information as function of radius implies that its mass is 
Mbh = (3.9 ±2.2) X 10^ Mq. The existence of intcrmcdiat c-mass black holes in globular 
clusters is consistent with several scenarios for globular cluster evolution proposed in 
the literature. The inferred mass for M15 is consistent with the extrapolation of the 
relation between Mbh and urms that has been established for galaxies. Therefore, 
these results may have important implications for our understanding of the evolution 
of globular clusters, the growth of black holes, the connection between globular cluster 
and galaxy formation, and the nature of the recently discovered 'ultra-luminous' X-ray 
sources in nearby galaxies. Instead of a single intermediate-mass black hole, M15 
could have a central concentration of dark remnants (e.g., neutron stars) due to mass 
segregation. However, we argue that the best-fitting Fokker-Planck models that have 
previously been constructed for M15 do not predict a central mass concentration that 
is sufficient to explain the observed kinematics. To fit the M15 data without any 
central dark mass concentration one must assume that the velocity distribution is 
significantly radially anisotropic near the center, which contradicts predictions from 
both Fokker-Planck and AT-body calculations. 

Subject headings: globular clusters: individual (M15) — stars: kinematics. 
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1. Introduction 

The globular cluster M15 (NGC 7078) has one of the highest central densities of any globular 
cluster in our Galaxy. As a result, it has been one of the globular clusters for which the structure 
and dynamics have been most intensively studied in the past decade (as reviewed in van der Marel 
2001). The present paper is the second in a series of two in which we present the results of a study 
with the Hubble Space Telescope (HST) of the line-of-sight velocities of stars in the central few 
arcsec of M15. Paper I (van der Marel et al. 2002) discussed the observations, and the extraction 
and calibration of the stellar spectra. In the present paper we determine the stellar line-of-sight 
velocities from the spectra, and we use the results to study the dynamics and structure of M15. 

M15 is a proto-typical core-collapsed cluster (Djorgovski & King 1986; Lugger et al. 1987; 
Tragcr, King & Djorgovski 1995), with a stellar surface density profile that rises all the way 
into the center. Such clusters make up ~ 20% of all globular clusters in our Galaxy, and stand 
in marked contrast to King-model clusters, which show fiat central cores and are modeled as 
tidally-truncated isothermal systems. Even imaging studies with HST have not provided any 
evidence for a homogeneous core in M15 (despite early claims to the contrary; Laucr et al. 1991). 
Guhathakurta et al. (1996) used the Second Wide Field and Planetary Camera (WFPC2) and 
found the projected surface number density profile inside 6" (0.34 pc) to be consistent with a 
power law N{R) oc i?-0-82±o.i2 g^g-^ ^ j^j^^g ^^ggj^ ^ggjj ^Yie Faint Object Camera (FOG) and 
obtained N{R) oc iJ-O-'i'Oio.os £qj, turnoff stars. They also showed that the distributions for stars 
of different masses have slightly different power-law slopes, which is qualitatively consistent with 
the mass segregation predicted in a cluster in which two-body relaxation has been important. 

Bahcall &: Wolf (1976, 1977) constructed detailed models for the equilibrium stellar density 
distribution of a globular cluster in which a central black hole (BH) has been present for 
much longer than the two-body relaxation time. For a cluster of equal-mass stars one expects 
N{R) oc in surprisingly good agreement with the observed star count profile for M15. 

While BHs have been convincingly detected in the centers of galaxies (e.g., Kormendy & Gebhardt 
2001), no convincing detections exist for globular clusters. On the other hand, few, if any, previous 
studies have had sufficient sensitivity to unambiguously detect BHs in globular clusters with 
masses Mbh ^ 5 x 10^ Mq. There are many ways in which globular cluster evolution at high 
densities can lead to the formation of a massive BH in the center (Rees 1984). For example, core 
collapse induced by two-body relaxation may lead to sufficiently high densities for individual stars 
or stellar-mass black holes to interact or collide, with a single massive BH as the likely end product 
(Sanders 1970; Quinlan & Shapiro 1987, 1990; Lee 1987, 1993, 1995). Studies of such scenarios 
have gained much interest lately (Miller & Hamilton 2002; Mouri & Taniguchi 2002; Portegies 
Zwart & McMillan 2002) after the discovery of intermediate luminosity X-ray objects in external 
galaxies (e.g., Colbert & Mushotzky 1999). The emission of these objects may be due to accretion 
onto intermediate mass BHs. However, this interpretation is not uniquely implied by the data and 
there is no unique association of these objects with star clusters (e.g., Zezas & Fabbiano 2002). 

While the observed star count profile of M15 is consistent with the presence of a BH, it can 
be explained equally well as a result of core-collapse (Grabhorn et al. 1992). Hence, the star count 
profile by itself yields only limited insight. An additional problem is that photometric studies 
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cannot determine whether Hght follows mass, and what the abundance and distribution of dark 
remnants are. Kinematical studies arc therefore essential to gain further insight. Integrated light 
measurements of M15 initially suggested a very high central dispersion, a = 25 ± 7kms~^. This 
was a sharp increase from the dispersions of 5-15 kms~^ found at larger radii from the radial 
velocities of individual stars, which was interpreted as evidence for the presence of a 10^ Mq 
central BH (Peterson, Seitzer & Cudworth 1989). This pioneering work spurred a lot of interest 
in globular cluster dynamics. However, all subsequent studies were unable to confirm the high 
central velocity dispersion measurement. In particular, it became clear that that the weighting of 
stars by their brightness in an integrated-light spectrum produces a large systematic uncertainty 
in the velocity dispersion deduced from the broadening of lines (Zaggia, Cappaccioli, & Piotto 
1993; Dubath et al. 1994). Velocity measurements of individual stars are therefore called for. 

Line-of-sight velocities are now known from ground-based studies for ~ 1800 M15 stars, 
as compiled by Gebhardt et al. (2000a). Many different studies contributed to this dataset, 
as reviewed in the introduction of Paper I. The projected velocity dispersion profile increases 
monotonically inwards from cr = 3 it 1 kms~^ at i? = 7 arcmin (Drukier et al. 1998), to 
£7 = 11 ± 1 kms~^ at R = 24!' . The analysis of Gebhardt ct al. (2000a) suggested that the velocity 
dispersion is approximately constant at smaller radii, and is o" = 11.7 it 2.8 km s~^ at the innermost 
available radius i? f» 1". Outside of the very center, the velocity dispersion profile is well fitted by 
a spherical dynamical model with an isotropic velocity distribution and a constant mass-to-light 
ratio T = 1.7 (in solar V-band units). However, this model underpredicts the velocity dispersion 
in the central 2 arcsec. The fit can be improved by addition of a central black hole, which causes 
the velocity dispersion to rise in Keplerian fashion as o" cx towards the center of the cluster. 

The best fit was obtained with a mass Mbh ~ 2 X 10^ Mq (Gebhardt et al. 2000a). However, the 
ground-based M15 velocity dispersion data can be fitted equally well with a model in which the 
mass-to- light ratio T(r) of the stellar population increases inwards to a value of ~ 3 in the center. 
This would not a priori be implausible, since mass segregation would tend to concentrate heavy 
dark remnants to the center of the cluster. Models with an anisotropic velocity distribution may 
even be able to fit the data with constant mass-to-light ratio and without a central black hole. 
Higher spatial resolution data and more detailed modeling are necessary to decide amongst these 
scenarios; this is the focus of the present series of papers. 

It has been known for some time (Gebhardt et al. 1994) that M15 has a net projected rotation 
amplitude of V^ot ^ 2 kms~^ at radii comparable to the half-light radius (about 1 arcmin). More 
recent work (Gebhardt et al. 2000a; these results were also suggested by the integrated-light 
measurements of Peterson 1993) has revealed that the rotation amplitude is larger at small radii: 
V-rot = 10.4 lb 2.7 km s"-*^ for R < 3.4", implying that K-ot/c ~ 1 in this region. This large amplitude 
is surprising because two-body relaxation should rapidly transfer net angular momentum outward 
from such small radii (see the discussion in Gebhardt et al. 2000a). Even more surprising is that 
the position angle of the projected rotation axis at small radii is ~ 100° different from that near 
the half-light radius. Although the large increase in the rotation amplitude at small radii may 
have something to do with the presence of a central BH (Gebhardt et al. 2000a), the increase and 
change in position angle are not predicted by any current theory of globular cluster dynamical 
evolution. 
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Phinney (1993) used an alternative argument to constrain the mass distribution of M15. 
There are two mihisecond pulsars in M15 at a distance R = 1.1" from the cluster center that have 
a negative period derivative P. This must be due to acceleration by the mean gravitational field 
of the cluster, since the pulsars axe expected to be spinning down intrinsically (positive P). The 
observed P values place a strict lower limit on the mass enclosed within a projected radius of 
R = 1.1" . Combined with the observed light profile this implies that the mass-to- light ratio must 
increase centrally inwards. A similar pulsar acceleration study was recently performed by d'Amico 
et al. (2002) for the cluster NGC 6752, which suggests a central increase in mass-to-light ratio in 
this cluster as well. For M15, Phinney (1993) obtained T > 2.1 for the total mass-to-light ratio 
within R < 1.1" , with a statistically most likely value of T ~ 3.0. These results arc consistent 
with the analysis of stellar kinematics (see also Dull et al. 1997). Unfortunately, the pulsar data, 
like the ground-based observations of the kinematics, does not constrain the distribution of mass 
tightly enough to discriminate between the effects of mass segregation and a central BH. 

Tighter constraints on the distribution of mass near the center of M15 need observations 
of the kinematics with better angular resolution than previous studies and observations from 
space can supply these. So we started a project to use HST to determine more stellar velocities 
close to the center of M15 (HST program GO-8262, PI: van der Marel). As described in Paper 
I, we used the Space Telescope Imaging Spectrograph (STIS) to obtain observations with the 
0.1"-wide slit at 18 adjacent positions near the cluster center. All spectra cover the wavelength 
range from 5073-5359A, which includes the Mg b triplet at ~ 5175A. The resolution is 0.276A per 
pixel, which corresponds to 15.86 km s~^. Extensive reductions and calibrations were performed 
to extract spectra with signal-to-noise ratio S/N > 5.5 per pixel for a total of 131 stars. The 
velocity calibration of the spectra was the most crucial and difficult aspect of the data reduction. 
Corrections were necessary for: drifts in the wavelength scale during an orbit; changes in the 
velocity of HST as it orbits the Earth; and wavelength shifts induced by the offsets of stars from 
the center of the slit. The analyses in Paper I indicate that the uncertainty in velocity scale caused 
by residual calibration errors in the final spectra is ~ 2.5 knis^-*^. 

Here we analyze the 131 stellar spectra from Paper I, and we show that for 64 of them the 
quality is sufficient to obtain an accurate line-of-sight velocity measurement. We use the results to 
obtain new constraints on the dynamics and structure of M15. The paper is organized as follows. 
In Section 2 we describe the cross-correlation algorithm that we have used for the extraction 
of line-of-sight velocities, including the choice of spectral templates. We discuss the reliability 
of the results based on an analysis of our STIS observations of a calibration star. In Section 3 
we describe the application of the cross-correlation algorithm to the STIS spectra of M15. We 
describe how we have corrected the inferred velocities for the effects of crowding and blending. 
The reliability of the inferred velocities is verified by comparison to ground-based data, for those 
stars for which the latter are available. In Section 4 we infer the velocity dispersion and rotation 
velocity profiles of M15 from the combined HST and ground-based line-of-sight velocity samples. 
In Section 5 we present dynamical models to interpret the results, and we discuss the implications 
for the dynamical structure and mass distribution of M15. Section 6 discusses and summarizes the 
main conclusions. 
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2. Extraction of Line-of-Sight Velocities 
2.1. Algorithm 

Extracting line-of-sight velocities from stellar absorption line data is usually accomplished 
by cross-correlating the spectral data with a template spectrum of known velocity. One of the 
most widely used implementations of the cross-correlation method has been developed by Tonry 
and Davis (1979). Several non cross-correlation based methods have also been developed in the 
past two decades. While these methods are better suited to extract the full line-of-sight velocity 
distributions from absorption line spectra, they do so at the cost of requiring high signal-to-noise 
ratio spectra. For the STIS data analyzed here, the only kinematical quantity of interest is the 
line-of-sight velocity. These data are thus best suited to a cross-correlation based analysis. 

To this end, all line-of-sight velocities were derived in IRAF using the task XCSAO. XCSAO is 
part of the RVSAO package described extensively in Kurtz &; Mink (1998). XCSAO is essentially 
a refined and updated version of the cross-correlation algorithm developed by Tonry & Davis 
(1979). XCSAO is a well-tested and often used task with many attractive features. It provides a 
handle on both the confidence and the error of a particular velocity measurement by means of the 
r statistic, developed originally by Tonry Sz Davis (1979), which we refer to here as Vcc (in essence, 
this statistic measures the amplitude of the cross-correlation peak, divided by the amplitude of 
the 'average' peak expected from noise and template mismatch). Following Kurtz & Mink (1998), 
we chose to calibrate this statistic empirically (see below). Another useful feature of XCSAO is its 
batch mode option which considerably simplifies the task of having to analyze large numbers of 
spectra (our STIS data set consists of 19200 apertures). 

All spectra were reduced and wavelength calibrated with an adaptation of the STIS reduction 
pipeline as described in Paper I. Before the actual cross-correlation takes place, XCSAO performs 
the following tasks: continuum removal, apodization and Fourier filtering. Due to undersampling 
of the line spread function (see Paper I) the STIS spectra are strongly undulated, i.e. there are 
low frequency flux variations with wavelength that are not related to the intrinsic continuum 
distribution. To roughly preserve the correct line ratios, the continuum was removed by division 
rather than by subtraction. The ends of the spectra were cosine tapered to avoid aliasing. Finally, 
a Fourier bandpass filter was applied to remove both the low frequencies (residual continuum 
variations) and the high frequencies (noise). The software then calculates the cross-correlation 
function and finds the strongest peak in a 200 kms~^ range centered on the systemic velocity of 
M15 (—107.5 kms~^). The position of this peak is fitted with a parabolic function to derive the 
velocity and its formal random error. 



2.2. Spectral Templates 

The strength of a cross-correlation peak depends largely on how well the template spectrum 
matches the observed spectrum. Template spectra were therefore obtained of a cluster star at 40 
arcsec from the center of M15, which was also used for target acquisition purposes (see Paper I). 
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This st ar is a very bright giant and relatively high signal-to-noise spectra were therefore readily 
obtained. This ensures that uncertainties in the inferred velocities result mainly from the noise in 
the M15 data and are not associated with the template (save template mismatch). 

Templates observed with the same instrumental setup as the science observations have the 
advantage that the instrumental signature will be similar in both spectra. Any instrumental 
effect can therefore be easily calibrated out. However, the largest instrumental effect is usually 
a broadening of the cross-correlation function due to instrumental broadening of the absorption 
lines. This does not influence the position of the cross-correlation peak and, hence, does not affect 
velocity measurements. 

In addition to the observed templates, a number of artificial templates were created based 
on Kurucz models. A grid of 30 model stellar atmospheres with varying temperatures (4000K to 
8000K), metallicities (—2.1 dex to —2.5 dex) and surface gravities (log 3 from 0.5 to 2.0) covering 
a wavelength range from 5000 to 5500 A were calculated. The metallicity range was chosen to 
bracket the current best estimate of M15's metallicity. 

Artificial templates have the benefit of having infinite signal-to-noise ratios and they can be 
constructed to closely match the spectral features of the M15 data. They lack the instrumental 
signature, but as noted above, this would mainly affect measurements of the velocity dispersion 
(i.e., the line widths) which is irrelevant to the data presented here. 



2.3. Reliability of the Line-of-Sight Velocities 

In order to successfully attain our goal of measuring the stellar velocity dispersion within 
the central arcseconds of M15, several intricate corrections to the derived line-of-sight velocities 
had to be made. These corrections stem from the motion of HST (which has an orbital velocity 
of 7.5 kms^^) and the fact that the stars generally do not fall in the center of the slit (the full 
width of the slit corresponds to 26.5 km s"^ in the dispersion direction). To calibrate the necessary 
corrections, 14 short exposure spectra of the bright field star HD 122563 (F8IV, magnitude 
V = 6.2) were obtained at offset positions parallel and perpendicular to the slit. Both the 
corrections and the calibration are described fully in Paper I. In the remainder of this paper it 
should be implicitly understood that these corrections were applied to all velocities. The star HD 
122563 was chosen because its low metallicity, [Fe/H] = —2.65 it 0.2 (Sneden &; Parthasarathy 
1983), which is close to the metallicity of M15, [Fe/H] -2.22 (Harris 1996). 

An additional use of the two-dimensional long-slit data obtained for HD 122563 is to calibrate 

the confidence limits on the results obtained with the XCSAO cross correlation software. In each 
of the 14 long-slit spectra we extracted 21 one-dimensional spectra at different positions along 
the slit. This yielded 294 spectra, each for a 0.1" x 0.1" aperture placed somewhere within 1.0" 
from the position of the star on the sky. Because of the wings of the PSF, each of these spectra 
contains a meaningful amount of light from the star. For each of these 294 spectra the line-of-sight 
velocities and the errors were derived using the cross-correlation routine XCSAO. In addition, the 
Tec statistic was derived, as well as the average signal-to-noise ratio S/N per pixel (determined 
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from the error frames delivered by the HST/STIS pipehnc). All apertures were correlated against 
all Kurucz model templates and also against the acquisition star template. The results were 
found not to change significantly from template to template. This is most likely due to the rather 
limited wavelength range of the observed spectra, 250 A. In the end we adopted the template 
which yielded the marginally highest value of Tcc- 

The inferred velocities were used to empirically calibrate the confidence limits in a manner 
analogous to Kurtz & Mink (1998). The left panel of Figure 1 shows the logarithm of the absolute 
difference between the measured line-of-sight velocity and the literature value of the velocity of 
HD 122563 versus the Vqc statistic. A clear break in the distribution of points occurs at v^c ~ 2.5. 
For Tec ^ 2.5 all velocity differences are only a few kms~^ and are thus considered reliable. Not 
surprisingly, for small values of the inferred velocities become unreliable. As an additional test 
of the reliability we have constructed a similar diagram showing the logarithm of the absolute 
velocity difference versus the signal-to- noise ratio (right panel of Figure 1). This plot shows the 
same qualitative behavior as the Vcc statistic plot. It suggests that reliable velocities are derived 
for an average S/N ;^ 5 per pixel. 

An optimal selection of the Tcc and S/N cutoff limits can be made by applying both criteria 
simultaneously. Imposing a lower limit on Vcc of 2.0 (vertical line in the left panel) would include 
very deviant points. However, all of the deviant points have S/N values smaller than 5.5 (these 
points are shown as filled circles in the left panel of Figure 1). A similar situation is encountered 
in the right panel, where all the deviant points that are beyond the S/N cutoff, S/N ^ 5.5, have 
unacceptable Tcc values (again shown as filled circles). Thus, in both panels all points that are to 
the right of their adopted cutoffs and are shown as open circles are considered acceptable. With 
the adopted values of = 2.0 and S/N = 5.5, no deviant points remain while the number of 
reliably measured velocities is maximized. Extending the and S/N cutoffs downward as far as 
possible is by no means critical for the calibration star itself but in M15 the number of reliably 
measured stars is at a premium. Kurtz &; Mink (1998) generally use more conservative values of 
Tec- However, they do not apply the Vcc statistic and the S/N criterion simultaneously to their 
data. 

The additional checks discussed below also indicate that the values adopted here yield reliable 
results for our STIS setup. 



3. Kinematical Results 
3.1. Raw Stellar Velocities 

As described in Paper I, the STIS data set allowed the extraction of (one-dimensional) spectra 

with S/N > 5.5 for a total of 131 stars in M15. We analyzed each of the extracted one-dimensional 
spectra with the cross-correlation algorithm described in Section 2. Based on the analysis of the 
calibration star discussed in Section 2, we expect that the inferred velocities are reliable for all 
spectra that yield a cross-correlation statistic rcc > 2. As a sanity check, we plot in Figure 2 the 
quantity {v - Vgys) /cr as a function of Tcq. Here fsys 

—107.5 kms ^ is the systemic heliocentric 
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velocity of M15 and cr ~ 12 kms~ is a rough estimate of the stellar velocity dispersion near the 
center (Gebhardt et al. 2000a). This figure is analogous to Figure 1. It shows that there are no 
inferred velocities with > 2 that are implausible in view of our understanding of the kinematics 
of M15. So there seems no reason to mistrust any of the inferred velocities with Tcc > 2, which 
confirms the results from the calibration star analysis. 

The STIS data yield a reliable velocity result (i.e., S/N > 5.5 and rcc > 2) for 64 of the 131 
extracted one-dimensional spectra. According to the prescription of Paper I, the 131 extracted 
spectra were constructed by combining one or more apertures. It makes essentially no difference, 
however, whether the velocities are derived directly from these co-added spectra, or whether 
they are derived by averaging the velocities of individually analyzed apertures. The 64 selected 
stars are listed in Table 1. We will refer to this sample as our 'STIS velocity sample'. The 
identification number in the first column corresponds to the entry number in the HST/WFPC2 
stellar catalog described and presented in Paper I. The inferred heliocentric velocity Wobs is listed 
in column 6, and its error Ai^obs in column 7. These are the 'raw' velocities inferred from the 
spectra, which are not yet corrected for blending. Blending corrections are discussed in Section 
3.3 below. The velocity error Auobs listed for each star is the quadratic sum of the formal error 
of the cross-correlation result and the spectral velocity calibration uncertainty. For the latter we 
adopted a value of 2.5 kms~^, based on the results and discussion of Paper I. The analysis of the 
calibration star data shown in Figure 1 confirms that there there are no systematic uncertainties 
in the data that exceed this value. 

Figure 3 shows the color- magnitude diagram (CMD) of M15 for the stars in the HST/WFPC2 
catalog discussed in Paper I. The stars in the STIS velocity sample are highlighted. These stars 
generally have a magnitude in V brighter than 18.5, although there is one star that is as faint as 
V = 19.05. Red giant and sub-giant stars dominate the sample. This is true not only because their 
number density is higher than for bluer horizontal branch stars, but also because blue stars often 
lack sufficiently strong stellar absorption lines to yield a strong cross-correlation peak (Mayor 
1980). Two stars in the STIS velocity sample (ID numbers 3393 and 5846; encircled crosses in 
Figure 3) occupy a location close to the instability strip and are probably RR Lyrae variables. RR 
Lyrae stars are velocity variables that can have velocity excursions of zb50kms~^ from the mean 
(Smith 1995). We therefore exclude these two stars from our subsequent dynamical analysis. 

The distance d of the stars in the STIS velocity sample from the cluster center ranges from 
0.2 to 27 arcsec. The adopted position of the cluster center is given in Paper I, and is based 
on the analysis of Guhathakurta et al. (1996). The influence of the uncertainty in this position 
(0.2" in each coordinate) on the results of our study is negligible, as discussed in Section 5.3. 
The solid line in the left panel of Figure 4 shows a histogram of the distribution of the stars in 
the STIS velocity sample with distance. For comparison, the dashed line shows the distribution 
of the stars for which a velocity determination (with an uncertainty better than 10 kms~^) was 
known from ground-based data (Gebhardt et al. 2000a). The logarithm of their corresponding 
cumulative distributions are shown in the right panel. The heavy solid line in this panel shows the 
distribution for the combined STIS and ground-based velocity samples. It is clear that the STIS 
data significantly increase the number of stars with known velocities close to the cluster center. 
Figure 5 shows the spatial distribution in the central 4" x 4" of the stars with known line-of-sight 
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velocities. 



3.2. Consistency Checks 

About one third of tiie stars in the STIS velocity sample also have a ground-based measurement 
of the line-of-sight velocity (Gebhardt et al. 2000a). A comparison of the ground-based velocities 
to the HST/STIS velocities is shown in Figure 6. The agreement is generally excellent, which 
confirms the reliability of the HST results. The residual velocities (whst — ^'ground) are shown in 
Figure 7 as a function of F-band magnitude and B — V color. There are no trends with either, 
which once more validates the velocities determined with STIS. The weighted mean residual 
velocity is 0.2 it 1.4km s~^, indicating that our velocity zero point is consistent with what has been 
inferred from ground-based observations. The value is 29.1 for 23 points, which is statistically 
acceptable. 

Neither the HST velocities nor the ground-based velocities in Figure 6 contain an explicit 
correction for blending^. However, even the best ground-based data (Gebhardt et al. 2000a) has 

a poorer PSF than the HST data, so some of the small discrepancies seen in Figure 6 may be 
due to the effects of blending. This is true in particular for the two labeled stars in Figure 6, 
which appear to lie somewhat offset from the one-to-one relation. The first star, ID 5222, is close 
(~ 0.1") to another star and it is possible that the ground-based determination of its velocity is 
more strongly affected by blending than the HST measurement. The second star, ID 5831, is a 
member of the triple star system AC 214 that is located within 0.5" of the M15 center. Its three 
V 15.5 members are separated by about 0.1". Our STIS sample contains the velocities of two 
of AC 214's members (ID 5831 and ID 5846). We were unable to extract a reliable velocity for 
the bluest of the three stars (ID 5872, B-V = 0.02). Interestingly, for ID 5831 and ID 5846 we 
find a velocity difference of ~ 50 kms~^. This provides a tantalizing hint for the existence of a 
central compact object in the midst of AC 214. However, ID 5846 is a possible RR Lyrae variable 
and the measured velocity difference is consistent with the maximum velocity excursion of an RR 
Lyrae variable. The ground-based measurement (— 116kms~^, Gebhardt et al. 2000a) of ID 5831 
is consistent with being the weighted mean of the two values inferred from the HST/STIS data. 
The ground-based measurement therefore probably suffered from blending of light with ID 5846. 

In addition to the comparison with ground-based results, we have also performed internal 
consistency checks on the STIS data set. As described in Paper I, the observational strategy 
yielded several independent spectra for most stars. These independent spectra correspond either 
to different rows in the same exposure (i.e., adjacent positions along the slit for a given slit 
position) or to different exposures (i.e., adjacent slit positions on the sky, observed either within 
the same telescope visit to the target or in different visits). The strategy adopted for our final 
results has been to co-add spectra for a given star, and to analyze the grand-total summed 



^The ground-based velocities cannot easily be corrected for blending. So when comparing the ground-based 
velocities to the HST velocities it makes sense to use HST velocities that are not corrected for blending either. In 
the subsequent dynamical analysis we do correct the HST velocities for blending, as described in Section 3.3 below. 
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spectrum. However, for testing we have also performed a cross-correlation analysis of all 19.200 
individual spectra in our dataset. We extracted those results that conform to the selection criteria 
with respect to blending, S/N and Tec that were discussed previously. For many stars in our STIS 
velocity sample this yields two or more independent velocity measurements. Let Vi and Vj be two 
independent velocity measurements for the same star, and let Avi and Avj be the corresponding 
uncertainties. One then expects each difference Avij = Vi — vj to be a Gaussian random deviate 
drawn from a distribution with zero mean and dispersion eij = [(Afj)^ + (Auj)^)]^/^. Figure 8 
shows a histogram of all the Avij. It includes those stars for which multiple independent 
measurements could be extracted; each star can appear many times in the histogram, depending 
on how many different measurements are available. The histogram is well fit by a Gaussian with a 
dispersion of 6kms^^ (solid curve). This is in reasonable agreement with the expected value of 
5 kms~^, based of the average of the individual uncertainties e^j. This implies that there is no 
reason to suspect the presence of unidentified systematic errors in the analysis, consistent with the 
results from the comparison of STIS and ground-based results. 



3.3. Statistical Corrections for Blending 

Even with the resolution of HST, the light profiles of stars in the densely populated central 
region of M15 overlap partially with neighboring stars. So when the spectrum of an individual star 
is extracted from our HST/STIS data set, only a fraction / of the total light is actually due to the 
star of interest. The remaining fraction (1 — /) is due to 'contaminating' stars. The algorithm for 
spectral extraction employed in Paper I was designed to yield values / > 0.75. We verified that 
the main results of our analysis are insensitive to the exact choice of this cut (see also Section 4 
below). The fraction / is listed in Table 1 for each spectrum that was analyzed. The blending 
fraction (1 — /) < 0.25 is small, but does need to be corrected for. 

To correct for blending we note that the observed velocity Vohs (the velocity extracted from 
the one-dimensional spectrum using cross-correlation) is, to a good approximation, equal to the 
luminosity weighted mean velocity of all the light that contributes to the spectrum. Therefore, 

^^obs = fVstur + (1 - /)^^cont, (1) 

where Vgtax the velocity of the star of interest, and t^cont is the luminosity weighted mean velocity 
of the contaminating stars. The analysis of Paper I tells us which stars are contaminating the 
spectrum and what fraction of the light they contribute. However, we do not generally know 
the velocities of the contaminating stars. On the other hand, we do know something about 
the contaminating stars: they belong to the cluster. Hence, the velocities of the contaminating 
stars are random deviates drawn from the line-of-sight velocity distribution of the cluster. The 
expectation value of the statistical quantity Vcont is therefore equal to Vgys- The dispersion is equal 
to Aa, where a is the line-of-sight velocity dispersion of the cluster at the position of interest. The 
constant A < 1 depends on the number of stars that provide the contaminating light. If there is 
only a single star, then A = 1. Upon taking the first and second moments of equation (1) we can 
calculate the expectation value of i^star, which we denote i;*, and the dispersion of ^^star, which we 
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denote Av*. This yields 

V* = Vsys + ybobs - ^^sys], = jii^Vo^.f + (1 " f^A^a^^^, (2) 

where Avobs is the formal error of the velocity inferred from the observed spectrum through 
cross-correlation. In the following we assume that A = 1; this is conservative in the sense that it 
yields the largest uncertainty Av^. 

The quantities v^. and Au* are the final estimate of the stellar velocity and its error. They 
result from correcting the measurement v^hs ^ind its error Af obs for blending. The accuracy of 
these corrections was verified through Monte-Carlo simulations. Equation (2) shows that the 
correction is such that it increases the deviation of the stellar velocity from the cluster mean by a 
factor 1//. This is because blending will, on average, tend to draw measured velocities towards 
the cluster mean. Equation (2) also shows that the error Av^ of the corrected velocity is larger 
than the error Av^hs of the measured velocity. This is so because of our poor knowledge of the 
velocities of the contaminating stars. The values of and Av* are listed in Table 1 for all the 
stars in the STIS velocity sample. They do not depend sensitively on the exact choices of Vsys and 
a, for which we adopted Vsys — -107.5 kms"^ and a = 12kms~^ (Gebhardt et al. 2000a). 

The median value of / in our velocity sample is 0.85. The uncorrected velocities have a median 
l^obs ~ Vgysl of 8.8 kms~^ and a median error Awobs of 3.8 kms~^. The corrected velocities have a 
l^'* — ^^sysl of 11.0 kms~^ and a median error Av^, of 5.4 kms~^. Thus, the blending corrections are 
always relatively small. This is because (1 — /) is always less than 25 percent; stars in our data 
set with larger amounts of blending simply did not make it into our final STIS velocity sample. 



4. Rotation Velocity and Velocity Dispersion Profiles 

The kinematical analysis of the data and the subsequent dynamical modeling benefit from 
having the largest possible sample. We have therefore combined the 64 stars in the STIS velocity 
sample with the 1777 stars in the ground-based sample used previously by Gebhardt et al. (2000a). 
There are 23 stars in common to both samples (see Table 1 and Figures 6 and 7). Those stars were 
assigned a velocity based on the weighted mean of the two measurements^. One could argue that 
the STIS velocities should be more accurate, so that it might be more appropriate to ignore the 
ground-based velocity measurements for these stars. However, since the HST and ground-based 
velocities are generally in good agreement (Figure 7), this doesn't affect any of our final results. 
Based on their locations in the CMD diagram, a total of 19 stars (including two in the HST 
sample) are possible RR Lyrae candidates. These stars were excluded from the sample, yielding a 



*The 1-a error in a weighted average of A'" measurements Zi is given by l/{Azi)'^]^^^'^ . However, this 

underestimates the true error if there are systematic errors in addition to Gaussian random errors. These cases 
manifest themselves by having a that is larger than the number of degrees of freedom (where the is defined to 
measures the residuals of the data with respect to their weighted average). In these cases we increased the errors to 
yield a x^ equal to the number of degrees of freedom. 
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final sample of 1797 stars.^ 

We have used two non-parametric techniques to estimate kinematical profiles from the discrete 
velocity measurements: LDWESS and SROTATE. These techniques yield the velocity dispersion and 
the amplitude and position angle of rotation, as functions of radius. LDWESS is a smoothing-based 
technique that uses a locally weighted regression algorithm, and is described fully in Gebhardt 
et al. (1994). It finds a smooth curve for the kinematical profiles as a function of radius. At 
each radius, it uses a fixed fraction of the data points in both radial directions. The data points 
are weighted depending on their distance from the radius of interest using an inverted parabolic 
function. The smoothing parameter is the fraction of the data included in the fit at any given 
radius, for which we generally use 10%. The kinematical parameters are estimated using a 
maximum likelihood estimator that includes both the uncertainties on the individual velocity 
measurements and their assigned weights. The statistical uncertainties on the inferred profiles are 
calculated using Monte Carlo simulation. An undesirable feature of LDWESS is that it tends to wash 
out small scale features, which can bias the kinematical profiles close to the center. This problem 
is circumvented by the technique SROTATE, which is described fully in Gebhardt et al. (2000a) (the 
technique is not referred to as SROTATE in that paper, but it is described in the fifth paragraph of 
its section 3.1). It is similar to LDWESS in most respects, and differs primarily in that it allows for 
a variable size of the data window (i.e., of the smoothing). At most radii, it used a fixed fraction 
of the data points in both radial directions, as does LOWESS. However, at the boundaries of the 
dataset the window size shrinks to a smaller number of data points, although no smaller than a 
preset 'minimum window size'. This approach preserves sharp radial variations close to the cluster 
center, while also maintaining the desirable feature of LDWESS of lessening large variations due to 
discrepant data points. The results described in the remainder of this section were all obtained 
with SROTATE and a minimum window size of 11 data points. The main features of the results 
were found not to depend sensitively on the exact choice of the minimum window size. As a check 
we also determined the kinematical profiles with LDWESS. Although the results from LDWESS are 
somewhat unreliable at the smallest radii, it did yield results that are statistically consistent with 
those derived from SROTATE. Overall, the kinematical profiles presented here are robust, and do 
not depend sensitively on either the choice of the analysis algorithm of the choice of the algorithm 
parameters. 

Figure 9a-c show the radial profiles of the velocity dispersion a, the rotation velocity Vrot, 
and the position angle PAkin of the kinematical major axis. The 68.3% confidence bands are 
indicated with dotted curves. At radii larger than a few arcsec the data set is dominated by 
stars with ground-based velocity measurements and, as a consequence, the results derived at 
these radii are almost identical to those obtained previously by Gebhardt et al. (2000a, their 
Figure 13). Therefore, the results at the smallest radii are of primary interest here. Gebhardt et 



^Wc note that the ground-based velocities in the sample were not corrected for blending in the same way that the 
HST/STIS velocities were. This should not be a problem, because the ground-based velocities dominate the sample 
only at radii where crowding is less of a problem. We verified this with a number of tests, for example, by studying 
the kinematics inferred from stars of different magnitudes. From these tests we conclude that potential blending in 
the ground-based data does not affect the main conclusions of our paper at a significant level. 
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al. (2000a) found an increase in the rotation velocity of M15 towards the center of the cluster, 
with a maximum of Vrot ~ 10 kms~^ at a radius R ~ 2.5". We find that the rotation velocity 
actually continues to rise inward to a value of T^ot ~ 13 kms~^ at our innermost (trustworthy) 
radius, R 0.5". The velocity dispersion a shows a mild inwards decrease at small radii, from 
(T 14 kms"-*^ at i? 4" to a 10 kms"-*^ at i? 0.5". This is the same radial range over which 
the rotation velocity shows a strong inwards increase. The value of PA^in appears to decrease by 
~ 50° when going inward from i? ~ 4" to i? ~ 0.5", but this decrease is only marginally significant. 
The fact that V^ot/c" ~ 1 near the center of M15 is quite surprising. This implications of this were 
already discussed at length by Gebhardt et al. (2000a), and we will therefore not comment further 
on it here. 

The pressure in a stellar system is determined by the mean squared velocity of the stars 
(Binney & Tremaine 1987). In hydrostatic equilibrium, the pressure gradient balances the 
gravitational force. The quantity that is of primary importance from a dynamical point of view is 
therefore neither Vrot nor a, but the mean squared projected line-of-sight velocity V^^t + We 
define o-rms{R) = (^ot + '^'^V^'^j where the angle brackets denote the average over a ring of radius 
R on the projected plane of the sky. The most convenient way to determine (trms from the discrete 
velocity measurements is to use SROTATE while keeping the rotation velocity fixed at zero. The 
resulting profile is shown in Figure 9d. Figure 10 is similar, but for visualization purposes it also 
shows the individually measured stellar velocities for comparison. The quantity cjrms measures 
the spread among these individual velocities (properly corrected for the small contribution from 
measurement uncertainties). The uncertainty in the determination of cjrms scales as the inverse 
square-root of the number of data points, and is at all radii much smaller than urms itself. Outside 
the central few arcsec, the data are dominated by stars with ground-based velocity measurements 
and the (Trms profile is almost identical to that inferred by Gebhardt et al. (2000a). However, at 
smaller radii our results are more accurate because of the addition of velocity measurements from 
HST. We find that ctrms increases from ~ lOkms^^ at i? = 10" to ^ 14kms~^ at R = 3". Inside 
of the central few arcsec ctrms is consistent with being constant at ~ 14 kms~^. This exceeds the 
central values of 10-12 kms"^ inferred by Dull et al. (1997) and Gebhardt et al. (2000a). 

In our analysis we have included only those STIS spectra for which a single star contributes 
a fraction / of the total light that exceeds /min = 0.75. If one chooses a more stringent limit, 
then the number of stars for which velocities can be determined decreases. For example, if one 
uses /min = 0.90, then the size of the STIS velocity sample reduces to only 22 stars (see Table 1). 
On the other hand, these stars are the ones for which the measurements are most secure. It is 
therefore a useful test to check that the results for these stars are consistent with the results for 
the sample as a whole. The STIS velocity sample contains 13 stars with / > 0.90 that reside 
at i? < 6" from the M15 center. These stars have an average distance from the cluster center 
(R) = 3.3" and their RMS projected line-of-sight velocity is ctrms = 14.1 it 3.2 kms^"*^. This 
is in excellent agreement with the results obtained from the full sample, shown in Figure 9d. 
This indicates that main results of our study do not depend significantly on either the amount 
of blending allowed during the extraction of the STIS spectra (see Paper I) or the subsequent 
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blending corrections (see Section 3.3).^ 

To assess the implications of our new observational results we proceed in Section 5 by 
constructing detailed dynamical models. As discussed below, these models use a maximum 
likelihood approach that operates directly on the individually observed stellar velocities. This has 
the advantage that there is no dependence on arbitrary binning or smoothing parameters, such as 
the data window size in SROTATE. Nonetheless, the kinematical profiles in Figure 9 are important 
for providing insight into both the dynamical state of M15 and the results of the likelihood analysis 
(although the profiles are not themselves used in the dynamical modeling). 



5. Dynamical Modeling 
5.1. Formalism 

To determine the mass distribution in M15 we assume that the cluster is spherical and that 
it is in hydrostatic equilibrium. The assumption of sphericity is adequate because the isodensity 
contours of M15 are not far from circular, at least at the small radii that are of primary interest 
here (Guhathakurta et al. 1996). The assumption of hydrostatic equilibrium is also adequate, 
despite the fact that a cluster like M15 undergoes secular evolution as a result of two-body 
relaxation. The two-body relaxation time significantly exceeds the dynamical time and the 
evolution therefore proceeds through a scries of states that are all approximately in hydrostatic 
equilibrium. A spherical system in hydrostatic equilibrium obeys the Jeans equation (Binney k. 
Tremaine 1987): 

p dr r dr 

Here p is the mass density, $ is the gravitational potential, /3 = 1 — erf /cr^ is a measure of the 
anisotropy of the stellar velocity distribution, and af. and are the mean-squared velocities in 
the radial and tangential directions, respectively. In general, all these quantities are a function 
of the three-dimensional radius r. Note that the assumption of sphericity does not imply that 
rotation is fully ignored. It is possible to construct spherical dynamical models that rotate. 
Rotation contributes to the mean-squared velocities in the azimuthal direction, and therefore is 
included in equation (3). So while we do not model the rotation of M15 explicitly, we do include 
its contribution to the hydrostatic support. 

At a given position in three-dimensional space, the mean-squared velocity in the direction 



The algorithm of Paper I that constructs a grand-total spectrum for each star has more parameters than just 
/min (although /min IS Certainly the most important one). For example, an individual aperture is rejected in the 
construction of the grand-total spectrum for a star if that star does not contribute at least a fraction g > pmin = 0.5 
of the light in the aperture. As for /min, the results of our analysis are not sensitive to the precise choice of gmin- The 
reason for this is that most of the light in the final grand-total spectrum for each star always comes from apertures 
with g > /min (otherwise it would not be possible for the sum of the aperture spectra to have / > /min). Since we 
choose /min S> gmin, the precisc choice of gtmin has negligible influence on the final spectra. 
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along the line-of-sight is af = a^[l — /3(r)i?^/r^], where R is the projected radius on the sky. 
The calculation of the observable projected kincmatical quantities depends on whether one 
observes integrated light or individual stars. In the integrated light situation, the mean-squared 
projected line-of-sight velocity (t|^ms line-of-sight projection of jaf, divided by the 

line-of-sight projection of j. Here j{r) is the three-dimensional luminosity density, which upon 
line-of-sight projection yields I{R), the projected intensity. By contrast, for M15 we are in 
the discrete situation, which is generally more complicated. Instead of j{r) one must know 
the three-dimensional number density of that population of stars from which the stars with 
kinematical measurements are drawn. This quantity is more difficult to estimate. To eliminate 
this complication we assume throughout our analysis that the number density of these stars 
is linearly proportional to the luminosity density j(r), with a proportionality constant that is 
independent of position. So we assume that the average luminosity per observed star does not 
depend on position. This is a minor simplification, because it doesn't allow for gradients in the 
visible stellar population within the cluster. Note, however, that no assumption is made about 
the distribution of dark objects in the cluster, which is the much more important quantity for the 
dynamical predictions of the modeling. 

In practice, we start the modeling from the V-band surface brightness profile ^y{R) in 
mag/arcsec^. For this we take the profile compiled from various sources by Trager et al. (1995). 
We correct this profile for Ay = 0.28 mag of extinction, which is based on the reddening 
E{B -V) = 0.09 (Harris 1996) and a standard Galactic reddening law with Av/E{B -V) = 3.1 
(Binney & Merrifield 1998). We transform this profile into a projected intensity I{R) in pc~^. 
The Trager et al. profile is based exclusively on ground-based data. Because of this we do not use 
their profile inside oi R = 5", where it is compromised by seeing and crowding. For R < 5" we use 
the projected number density profile derived non-parametrically by Guhathakurta et al. (1996) 
from HST/WFPC2 data, scaled to match the Trager et al profile at i? = 5". This profile rises into 
the center of M15 approximately as a power law, I{R) oc i?^*^'^^. Sosin & King (1997) inferred 
a slightly different profile from HST/FOC data, and they also found slightly different power-law 
slopes for different stellar masses. Inside ~ 0.3" the number density profile is not well constrained 
by either study, due to the limited number of stars and uncertainties in the position of the M15 
center. We have performed a variety of tests to assess the influence of these uncertainties on the 
predictions of our dynamical models. We have found that the main results are quite robust and 
do not depend sensitively on the exact choice of I{R) in the central 5". 

To calculate isotropic (/3 = 0) models for M15 we have used the non-parametric technique 
of Gebhardt &; Fischer (1995). It solves an Abel integral to calculate j{r) from I(R). The mass 
density is then calculated as p{r) = T(r)j(r), where T(r) is an assumed mass-to-light ratio profile. 
Poisson's equation is then solved to obtain the gravitational potential $(r). To this we add the 
contribution —GMbu/t of a possible central massive black hole. The Jeans equation (3) is then 
solved, and the results are projected along the line of sight to obtain (Trms(-R)- To calculate 
anisotropic models (/3 7^ 0), we have used the technique of van der Marel (1994). This technique 
starts with a parametric fit to I{R), but otherwise proceeds similarly as the technique of Gebhardt 
&: Fischer (1995). As a test, we used both techniques to calculate predictions for the case (3 = 
and T(r) = constant, which yielded excellent agreement. In all model calculations we assumed a 
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distance to M15 of 10.0 kpc (Harris 1996). 



5.2. Comparing the Data and Models 

Data-model comparison can be performed qualitatively by visual comparison of a predicted 
profile (Trms(-R) to the profile that was inferred from the data and is shown in Figure 9. However, 
quantitative data-model comparison is best performed using a maximum likelihood approach. 
We assume that the predicted velocity distribution at projected radius i? is a Gaussian with 
dispersion o"rms(^)- This is a reasonable assumption, at least for isotropic models (e.g., van der 
Marel & Franx 1993). Calculation of the full line-of-sight velocity distribution shape for spherical 
anisotropic models is not straightforward (e.g., Gerhard 1993; van der Marel et al. 2000) and is 
outside of the scope of the present paper (see Section 5.5 for a further discussion of anisotropic 
models). Let a star, number i, be observed at projected radius Ri. Let the uncertainty in its 
observed velocity Vi be Avi. The probability distribution Pi{v) for the velocity v observed for star 
i is then a Gaussian with dispersion [o"p^Ms(^«) + Avf]^^"^ , centered on the cluster systemic velocity 
Vgys- The total likelihood L of the dataset is the product of the quantities Pi{vi) for all of the 
N = 1797 stars in the combined velocity sample, so that 

N N 

X = -2\nL = J2 ln27r[4Ms(^) + ^^f] + - ^'sys)V[4Ms(^i) + ^^fl- (4) 

i=l i=l 

The models are generally characterized by a certain set of parameters. The parameters that yield 
the minimum A value, Amim are the ones that provide the best fit. 

Once a best-fit model has been identified, there are two further questions that need to be 
answered: (a) is the best-fitting model statistically acceptable; and (b) what are the confidence 
regions around the best-fitting model parameters? To address the first question, we note that 
the last term in equation (4) is a quantity. If each Vi is drawn from a Gaussian distribution 
with dispersion [crms(^«) ~^ which is what we assume here, then this term follows a 

probability distribution with A'^ degrees of freedom. For large N, this distribution has a mean N 
and dispersion ■\/2N. So for the best-fit model to be statistically acceptable at the la level, the 
value Amin should be within ±-\/2N from the expectation value 

N 

(A) =J2ln27r[aius{Ri) + ^vh+N, (5) 

i=l 

where (Trms(-K) is the profile predicted by the best-fit model. To address the second question, we 
use a well-known theorem of mathematical statistics (e.g., Stuart Sz Ord 1991; used also by Merritt 
Sz Saha 1993) which states that the likelihood-ratio statistic A — Amin tends to a statistic in the 
limit of large N, with the number of degrees-of- freedom equal to the number of free parameters 
that have not yet been varied and chosen so as to optimize the fit. This is a consequence of the 
central limit theorem. Hence, the likelihood-ratio statistic A — Amin reduces to the well-known 
Ax^ statistic (e.g.. Press et al. 1992). This allows straightforward calculation of confidence regions 
around the best-fitting model parameters. 
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5.3. Isotropic Models with Constant Mass-to-Light Ratio 

We started by calculating isotropic models with a (y-band) mass-to-light ratio T that 
is independent of radius. We allow for the possibility of a central black hole of mass Mbh- 
Figure 11a shows the contours of the likelihood quantity A in the (T, Mbh) plane. The solid 
contours show the 1, 2, 3, and 4a confidence regions on the best-fitting model parameters. 
The best-fitting mass-to-light ratio is almost independent of Mbh- This is because the black 
hole only influences the dynamical predictions close to the center, so that T is determined 
more-or-less uniquely by the data at large radii. By marginalizing over Mbh, we infer a best-fitting 
mass-to-light ratio T = 1.6 it 0.1. For understanding the dynamical structure of M15, the more 
interesting quantity is the black hole mass. By marginalizing over T we infer a best-fitting value 
Mbh = (3.2 ± 2.2) x 10^ M©. The best-fitting model has hkelihood A = 13405. Equation (5) 
predicts that A should be in the range 13419 it 60, so the best-fit model is statistically acceptable. 

Figure lib compares the predicted and observed profiles of ctrms for fixed T = 1.6 and 
different black hole masses. The model with Mbh = predicts a decrease in ctrms towards the 
center. This is a direct consequence of the fact that the luminosity and number density of M15 
increase steeply towards the center (e.g., Tremaine et al. 1994). The observed ctrms profile does 
not show a decrease towards the center, requiring the presence of a central black hole for the 
models to fit the data. 

As discussed in Section 1, there are independent constraints on the central mass distribution 
of M15 from observations of pulsar accelerations. The measured accelerations of two pulsars at 
R = l.l" from the cluster center yield a lower limit Mcyi,min(l-l") on the M15 mass Mcyi(l.l") 
enclosed within a cylinder of this radius. Equation (3.6) of Phinney (1993), combined with the 
pulsar data in his table 2, yields Afiniii(l-1") = 4.5 x 10^ M©. The most likely (expectation) 
value for the enclosed mass, (Mcyi(l.l")), depends on the phase-space distribution of the pulsars 
in the cluster. Phinney (1993) quotes that for typical distributions Mcyi^min/ (Mcyi) = 0.6-0.75. 
This yields (Mcyi(l.l")) ~ (6.0-7.5) x 10^ Mq. For the best fit model in Figure 11 we find that 
Mcyi(l.l") = 6.3 X 10^ M0 (of which 3.2 x 10^ Mq in the form of a black hole, and 3.1 x 10^ Mq in 
the form of stellar objects). This is in excellent agreement with the pulsar constraints. 

In our stellar dynamical analysis we have used the position of the cluster center given in 
Paper I, which is based on the analysis of Guhathakurta et al. (1996). This position has a 1-a 
uncertainty of 0.2" in each coordinate, which corresponds to a two-dimensional 1-a uncertainty of 
0.3". This is small compared to the relevant scales in the analysis. For example, the kinematical 
profiles in Figure 9 only start at R = 0.5" (this is the average radius of the innermost 11 data 
points, where 11 is the adopted data window size in SROTATE). More importantly. Figure lib 
shows that models without a black hole fail to fit the data over the entire region R ^6". Shifts in 
the assumed position of the M15 center of order ~ 0.3" cannot change this. To quantify this, we 
repeated the likelihood analysis with several different assumed positions for the M15 center. These 
positions were obtained by shifting the nominal center position by ~ 0.3" in various directions. 
The resulting changes in the best-fitting Mbh were found to be ^ 0.1 x 10^ Mq. This is much 
smaller than the formal error on Mbh (2.2 x 10^ M0) and is therefore negligible. The same applies 
to the other models that are discussed in the remainder of this section. Hence, none of the main 
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results of the present study are affected by the uncertainties in the exact position of the M15 
center. 



5.4. Isotropic Models with Varying Mass-to-Light Ratio 

The assumption of a mass-to-light ratio T that does not vary with radius is obviously a 
considerable simplification. M15 does not have a constant density core, which suggests that 
core-collapse has occurred as a result of two-body relaxation. A natural consequence of two-body 
relaxation is mass segregation. In an attempt to reach equipartition of energy, heavy stars and 
dark remnants sink to the center of the cluster, which causes a central increase in T(r). This 
boosts the predicted velocities close to the center, which may obviate the need to invoke a central 
black hole in the models (Peterson et al. 1989). 

Detailed Fokker-Planck calculations are required to determine whether mass-segregation 
resulting from realistic initial conditions can fit the kinematical data for M15. We have not 
calculated such models ourselves. Instead, we have relied on the results published by Dull et 
al. (1997), which are the most recently published Fokker-Planck models for M15. The Dull et 
al. models were constructed to fit kinematical data at radii i? ;> 4". It is therefore natural to ask 
how these models compare to the new data that is now available for R <^ 4" . Dull et al. do not 
show the kinematical predictions of their models inside R = 3". However, they do show the run of 
T(r) in their best-fit model (their figure 12). We measured this profile from their paper, and used 
it as input into our Jeans models. We adopted the profile for the situation that is intermediate 
between a collapsed core and an expanded core, which provides the steepest central increase in 
T. We scaled T(r) uniformly downward by a factor 0.86, so as to best match the kinematical 
predictions shown by Dull et al. at radii R^3". This scaling allows for differences in the adopted 
distance, foreground absorption, and surface brightness profile. The resulting mass-to-light ratio 
profile does not vary much at the radii R ^ 0.5" for which kinematical data is available, and is 
between 1.4 and 2.5 throughout. However, at smaller radii T increases to values ^ 10, due to a 
central cusp that is strongly dominated by non-luminous neutron stars. 

We constructed hydrostatic equilibrium models with the T(r) profile thus obtained from Dull 
et al.^ and with a range of black hole masses. Figure 12b compares the predictions for ctrms to 
the observed profile from Figure 9d. Figure 12a shows the likelihood quantity A as a fTinction 
of MbH) together with the 1 and 2a confidence levels. The best-fitting model has likelihood 
A = 13412. Equation (5) predicts that A should be in the range 13417 ±60, so the best-fit model is 
statistically acceptable. The best-fitting black hole mass is Mbh = (4.5 it 2.1) x 10^ Mq. So even 
with a realistic mass-to-light ratio profile T(r) inferred from Fokker-Planck models, the models 
still rcqTiirc a black hole to fit the data. The black hole mass is not very different from the value 
derived in Section 5.3 under the assumption of a constant T. The reason for this is that the 
enclosed mass inside a (three-dimensional) radius of ~ 0.5" is only ~ 1.0 x 10^ Mq larger for the 



'^Figure 12 of Dull et al. (1997) shows the mass-to-light ratio profile after projection along the line-of-sight, and 
not the intrinsic three-dimensional mass-to-light ratio profile. We properly accounted for this in our modeling. 
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Fokker-Planck models than it is for the constant T models. So even though in the former models 
there is a central concentration of dark remnants, the total mass in these remnants is insufficient 
to explain the observed kinematics. In fact, with the T(r) profile from Dull et al. (1997), the 
models require a larger BH mass than for a constant T (namely, [4.5 it 2.1] x 10^ Mq versus 
[3.2 lb 2.2] X 10^ M0). This is because T(r) actually decreases radially inwards between ~ 1 arcmin 
and ~ 3" (presumably caused by mass segregation concentrating the giants in the center compared 
to less massive stars). The best-fit model in Figure 12 has a mass Mcyi(l.l") = 9.1 x 10^ Mq 
within a cylinder of radius R = 1.1" (of which 4.5 x 10^ M.q in the form of a black hole, and 
4.6 X 10^ Mq in the form of stellar objects). This is somewhat larger than the expectation value 
implied by the pulsar constraints, (Mcyi(l.l")) ^ (6.0-7.5) x 10^ Mq (Section 5.3), but the residual 
is not statistically unacceptable. 

A large uncertainty in Fokker-Planck models comes from the assumptions used for the 
retention factor of neutron stars. Dull et al. (1997) assume that all neutron stars are retained. 
However, the distribution of pulsar kick velocities suggests that most of the neutron stars 
(assuming that all pulsars come from neutron stars) should have been ejected from the cluster due 
to the low escape velocity of M15. The best estimates for the retention factor are generally around 
a few percent and never go above 10% (Drukier 1996; Davies & Hansen 1998; Pfahl, Rappaport 
&: Podsiadlowski 2002). Even more extreme, Hansen & Phinney (1997) suggest less than 1% of 
single neutron stars would be retained. Binaries, however, provide a mechanism in which the 
neutron star retention factor may be increased, but, even when assuming a large binary fraction, 
the retention factor is still below 10% (Davies & Hansen 1998; Pfahl, Rappaport & Podsiadlowski 
2002). Thus, it appears that the contribution from neutron stars may not be very significant for 
the central mass distribution in globular clusters. Heavy white dwarfs (1-1.3 Mq) do not suffer 
from these large birth kicks and so can provide some contribution to the central mass, but these 
objects are not as centrally-concentrated (Dull et al. 1997). The best way to understand the effect 
of stellar remnants is to include evolutionary models with realistic assumptions for the retention 
factor. The important conclusion in the present context is that the Dull et al. models strongly 
overestimate the expected neutron star retention, yet still fail to fit the M15 data without a central 
black hole. 

The number of stars that end up as neutron stars or white dwarfs in Fokker-Planck models 
depends critically on the assumed initial mass function (IMF). The more stars there are at the 
high end of the IMF, the higher the predicted central velocity dispersion at the present epoch 
(e.g., Grabhorn et al. 1992). However, there are good observational constraints on the present-day 
local mass function (MF) of M15, at least below the main-sequence turn-off (~ 0.7 Mq for M15). 
The IMF can therefore not be treated as a completely free parameter. The MF obtained from 
HST photometry by De Marchi & Paresce (1995) for a field at R = 4.6 arcmin from the M15 
center has a power-law slope x ~ 1.3 at masses above ~ 0.25 Mq (where x is defined so that a 
Salpeter MF has x = 1.35). Dull et al. (1997) use a global IMF slope x = 0.9 at these masses, 
which yields a local present-day MF slope x 1.2 at i? = 4.6 arcmin (figure 8 of Grabhorn et 
al. 1992). This indicates that the Dull et al. IMF is approximately correct, and if anything, is 
slightly too shallow. So there is little room for an increased number of high mass stars in the IMF, 
which would be needed to bring the Dull et al. models in better agreement with the kinematical 
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data presented here. One would actually like the Fokker-Planck models to reproduce also the 
observational constraints on the radial dependence of the present-day MF (De Marchi & Paresce 
1995; Sosin &; King 1997). A preliminary report of attempts at this (Murphy et al. 1997) indicates 
that such models predict a present-day total number of ~ 7000 non-luminous 1.4 Mq objects in 
M15. This is 30% less than the number of such objects in the Dull et al. models. Hence, more 
accurate modeling of the IMF is not likely to improve the fit to the kinematical data presented 
here. 

The Fokker-Planck models of Dull et al. (1997) are the most recent ones to have been 
constructed for M15, but they are not the only ones. The Dull et al. models built on previous work 
by Grabhorn et al. (1992). Both papers used the same Fokker-Planck methodology and adopted 
the same IMF for their best fit. The main difference is that Dull et al. had access to more modern 
and extensive data to compare their models to. The models of Dull et al. therefore supersede the 
models of Grabhorn et al. Also, the discussion that was presented above with respect to the Dull 
et al. models applies equally to the models of Grabhorn et al. In another study, Phinney (1993) 
fitted an independent set of Fokker-Planck models to M15 data. These models were calculated by 
Bryan Murphy using an extension of the code described by Murphy, Cohn & Hut (1990). These 
models had several nice features, including explicit inclusion of stellar evolution and neutron-star 
escape. However, the models used a pure power-law IMF. This is less accurate than the IMF 
used by Grabhorn et al. and Dull et al. and is not consistent with our current understanding of 
globular cluster IMFs (e.g., Paresce &; De Marchi 2000). As mentioned above, a model cannot be 
considered to be acceptable unless it fits observational constraints on the present-day MF. Also, 
Phinney's model that comes closest to fitting the presently available kinematical data for M15 
(his model zl20) must be evolved for 21.4 Gyr. This considerably exceeds the accepted age of the 
Universe (e.g., Lineweaver 1999). 

The above arguments indicate that there is an urgent need for more detailed Fokker-Planck 
modeling of the structure and evolution of M15. The models should include as much as possible 

of the relevant physics, including neutron star escape. The most current datasets must be used to 
constrain the cluster density, kinematics and IMF. With such a study it will be possible to explore 
in much more detail than has been possible here whether the observed kinematics of M15, both 
the central (Jrms and its dependence on radius, can be explained without requiring a central black 
hole. 



5.5. Anisotropic Models 

To further explore the range of models that can reproduce the observed kinematical profiles 
for M15 we have also constructed models with an anisotropic velocity distribution. Models with 
radial velocity dispersion anisotropy predict a larger projected velocity dispersion near the center 
of a stellar system than models with an isotropic velocity dispersion (e.g., Binney & Mamon 1982). 
Radially anisotropic models for M15 may therefore be able to fit the data without invoking a 
central black hole. To test how much radial anisotropy is required for this, we constructed models 
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in which the anisotropy function (3 = 1 — a^/a^ is parameterized as 

P{r) = Po + {Poc-(3o)^^. (6) 

J.Z _|_ q2 

This parameterization yields a convenient semi-analytical expression for the solution of the Jeans 
equation (van der Marel 1994). It corresponds to a profile that varies from /3 = /3o at the center to 
P = Poo at large radii. The transition occurs at the scale radius a. 

The vahic of Poo has little influence on the dynamical predictions near the center of the chistcr, 
which is where isotropic models without a central black hole fail to fit the data (cf. Figure 11). 
Without loss of generality we therefore kept Poo fixed at zero in our study. We also assumed that 
the mass-to-light ratio T is constant as a function of radius, and we did not include a black hole 
(i.e., Mbh = 0). This leaves three free parameters, Po, a and T, which were varied to optimize 
the fit to the data. The best-fitting model was found to have Pq = 0.65, a = 18" and T = 1.6. 
The inferred mass-to-light ratio is similar to that inferred for isotropic models (Section 5.3). 
The best-fitting model has likelihood A = 13401. Equation (5) predicts that A should be in the 
range 13395 it 60, so the model is statistically acceptable. Figure 13b compares the predictions 
for the best-fitting anisotropic model (dotted curve) to the observed urms from Figure 9d. For 
comparison, we also show the predictions for models in which Pq was kept fixed at Pq = 0, 0.2, 
0.4 and 0.8, and in which only a and T were varied to optimize the fit. Figure 13a shows the 
likelihood quantity A for these fits as a function of Pq, together with the 1 and 2a confidence levels. 
The best fit and its la error are Pq = 0.65 ± 0.2. 

The models in Figure 13 show that it is possible to fit the M15 data without invoking a 
central mass concentration. However, this requires that a.,- = l.dat (i.e., P ~ 0.65) throughout 
much of the central region of M15. Such a high value is difficult to reconcile with our current 
understanding of the dynamical structure of globular clusters. Even elliptical galaxies, which are 
collisionless, do not have such large anisotropics (e.g., Kronawitter et al. 2000). By contrast, in 
a cluster like M15 considerable two-body relaxation must have occurred near the center, which 
tends to isotropizc the velocity distribution (Peterson et al. 1989). This has been demonstrated 
explicitly, both with anisotropic Fokker-Planck calculations in energy-angular momentum space 
and with A/^-body calculations (Takahashi 1996; Baumgardt, Heggie & Hut 2002). In the region 
that contains the central ~ 20% of the cluster mass, the anisotropy parameter P remains between 
and 0.1 at all times during the cluster evolution. Larger radial anisotropics can develop in the 
outskirts of the cluster, but we verified explicitly that this does not change the interpretation of 
the central kinematics of M15. 



6. Discussion and Conclusions 

We have obtained high spatial resolution spectroscopy of the central region of the globular 
cluster M15 with the STIS spectrograph on board HST. The observational setup, calibration 
and spectral extraction were discussed in Paper I. Here we have analyzed the spectra with a 
cross-correlation technique to determine the line-of-sight velocities of individual stars. Our final 
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STIS velocity sample contains 64 stars. Two-thirds of the stars in this sample have their velocity 
measured for the first time. Half of the stars reside within a projected radius R = 2.4" from 
the center of M15. The new data set triples the number of stars with measured velocities in the 
central R < 1" of M15 and doubles the number in the central R < 2". Our analysis includes the 
necessary (small) corrections for the effects of blending with neighboring stars. Detailed tests 
on a calibration star and comparison to ground-based M15 data demonstrate that our velocities 
are accurate and trustworthy. We combined the STIS results with existing ground-based data 
to obtain a total sample of 1797 stars in M15 with known line-of-sight velocities. We use the 
combined sample to determine the radial profiles of the most important projected kinematical 
quantities: the rotation velocity V^ot- the position angle of the kinematical major axis, PAkin; the 
velocity dispersion, a; and the RMS velocity averaged over rings on the projected plane of the sky, 
CRMS- Our results differ from earlier work only in the central few arcsec. In particular, we find 
that ctrms rises to ~ 14 kms~^ at the innermost radii. This is somewhat higher than the value of 
10-12 km s~^ inferred previously from ground-based data (Dull et al. 1997; Gebhardt et al. 2000a). 

To interpret the results we constructed dynamical models based on the Jeans equation for a 
spherical system. We compared the model predictions to the data using a maximum-likelihood 
technique to obtain the best-fitting model parameters and their confidence regions. If the velocity 
distribution is isotropic, then M15 must have a central concentration of non-luminous material. 
This could be due to an intermediate-mass black hole. If one were reluctant to invoke such an 
object, then one alternative may be that M15 has a central collection of dark remnants (neutron 
stars and/or stellar mass black holes). This arises naturally in a globular cluster due to the mass 
segregation that occurs as two-body relaxation drives the system to equipartition. However, 
we argued that the best-fitting Fokker-Planck models that have previously been constructed for 
M15 (Dull et al. 1997) do not predict a large enough concentration of dark remnants to fit the 
data. It remains to be seen whether alternative Fokker-Planck models can be constructed that 
generate a more massive concentration of dark remnants from plausible initial conditions. It 
is useful to note in this context that an important uncertainty in Fokker-Planck models comes 
from the assumptions used for the retention factor of neutron stars. Dull et al. assumed that all 
neutron stars are retained, in contrast with most recent work which predicts that only ^ 10% 
will be retained. So Dull et al. may actually have overestimated the central concentration of dark 
remnants in their models. Another alternative scenario is to assume that deviations from isotropy 
in the velocity distribution may be responsible for the observed kinematics. However, to fit the 
kinematical data without any mass concentration one must assume that the velocity distribution 
is radially anisotropic near the center, (3o = 0.65 it 0.2. This contradicts the predictions of both 
Fokker-Planck models and A^-body calculations, which suggest that the velocity distribution in 
the central region of a globular clusters remains close to isotropic at all times during the cluster 
evolution. 

In view of the results that we have presented, the presence of an intermediate mass black 
hole in M15 appears to be the most plausible explanation of the data. As noted in Section 1, 
there are several mechanisms by which such a black hole could plausibly have formed. For the 
best-fit black hole mass we adopt the average of the values which were inferred in Sections 5.3 
and 5.4 using a constant mass-to-light ratio and the Dull et al. (1997) mass-to-light ratio profile. 
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respectively. This yields: Mbh = (3.9 ± 2.2) x 10^ Mq. This mass is consistent with the constraints 
on the central mass distribntion of M15 implied by observations of pnlsar accelerations (Phinney 
1993). The black hole mass inferred for M15 matches remarkably well with the understanding 
that has been developed for the presence of black holes in the centers of galaxies. For these black 
holes, there is a strong correlation between the black hole mass and the velocity dispersion of 
the bulge component (Gcbhardt et al. 2000b; Ferrarese & Merritt 2000). Figure 14 shows the 
available data points and the best fit from the recent compilation of Tremaine et al. (2002). For 
M15, the luminosity weighted mean velocity dispersion within the half light radius (1.06 arcmin, 
Harris 1996) is 12.1 kms""^ (this quantity was defined and calculated similarly as in Gebhardt et 
al. 2000b). At this dispersion, the estimated black hole mass fits perfectly on the extrapolation 
of the relation established for galaxies (Figure 14). Interestingly, a study of the globular cluster 
Gl in the Andromeda galaxy, performed simultaneously with the present study, has also provided 
evidence for a central black hole (Gebhardt, Rich &; Ho 2002). Like M15, this globular cluster 
fits perfectly on the relation shown in Figure 14. This independent research strengthens the 
interpretation of the M15 data in terms of an intermediate-mass black hole. It has generally been 
believed that globular clusters and galaxies form and evolve quite differently, so it could be that it 
is a mere coincidence that they fall on the same Mbh — cr relationship. However, it could also have 
some deep physical significance. For example, it may point to a new link between galaxy formation 
and globular cluster formation. Or it may point to a link between the black holes in these systems. 
For example, it could be that the massive black holes in galaxies grew from seed black holes that 
arose in clusters. There may also be a link with the intermediate luminosity X-ray objects that 
are known to exist in external galaxies, and which have been argued to be intermediate mass black 
holes. These issues will need to be explored with future observational and theoretical studies. 

Despite the interesting evidence for the presence of an intermediate mass black hole in M15, 
some words of caution are justified. All of the dynamical models that have been constructed for 
M15 remain somewhat idealized. This is true both for the Jeans models presented here and for 
the Fokker-Planck models presented elsewhere. For example, the Jeans models assume exact 
hydrostatic equilibrium, which is generally expected to be a good assumption (see Section 5.1). 
Nonetheless, during periods of particularly rapid evolution in the cluster structure this assumption 
could yield results that are biased. Fokker-Planck models can address the cluster evolution directly 
and do not need to rely on the assumption of hydrostatic equilibrium. On the other hand, the 
results of Fokker-Planck models depend strongly on the processes of stellar evolution and binary 
heating, both of which are generally modeled only in rudimentary ways. More generally, globular 
clusters are complicated systems from a theoretical viewpoint, much more so than galaxies, and 
not all of the essential physics may yet have been fully understood. There is some evidence 
from observations that this may indeed be the case. For example, M15 rotates quite rapidly in 
the central regions (cf. Figure 9), and this is not naturally explained by any theoretical model 
(Gebhardt et al. 2000a). Also, the ctrms of M15 appears to have a small dip at intermediate 
radii (at R ~ 13"), which is not naturally explained by any of the models constructed here 
(Figures 11-13). Further studies to test the observational reality of these features would be 
valuable, as would further theoretical work to address their origin. 

We have assumed throughout our study, as has previous work on M15, that the observed 
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kinematics are characteristic of the cluster, and are not contaminated by possible orbital motion 
of stars in binary systems. There are several reasons that make this a reasonable assTimption 
(Hut et al. 1992). The binary fraction of globular clusters is believed to be only of order 10%. 
Also, most of the stars in the velocity sample are red giants. Their relatively large radii imply 
that any binaries must have large separations and orbital velocities ^ 25 kms"^. Inclination, 
phase and ellipticity effects imply that for an average binary at a random epoch only a fraction of 
the velocity amplitude will be observed along the line of sight. These issues conspire to make it 
extremely challenging to identify even a few binaries in globular clusters from large line-of-sight 
velocity studies, even with high quality multi-epoch data (Pryor et al. 1989). One can turn this 
around to argue that the average observed kinematics of large samples of stars should not be 
influenced significantly by any orbital motion in binaries. This assumption is supported by the 
fact that the mass distribution inferred here from stellar kinematics agrees with that inferred from 
pulsar studies. Nonetheless, it would be useful for future studies to attempt a detailed quantitative 
assessment of the potential contamination of the line-of-sight kinematics of globular clusters by 
binaries. 

In the future it may be possible to strengthen the observational constraints on the central 
structure of M15 through proper motion studies. With two additional velocity components it will 
be possible to directly establish the (an)isotropy of the stellar velocity distribution. It has been 
demonstrated that such studies are feasible with HST (e.g., Anderson &; King 2000). However, 
the severe crowding in the central few arcsec of M15 may provide a significant hurdle to overcome. 
An alternative way to strengthen the observational constraints would be to increase the sample of 
radial velocities. This would require inclusion of fainter stars near the turnoff magnitude (V ~ 19 
in M15), which were inaccessible to our STIS study because of limited S/N. In principle, a 4m 
class telescope has sufficient light gathering power to perform such a study in a reasonable amount 
of time. By centering the spectra around the CO band-head (2.3/im), it may be possible to take 
full advantage of adaptive optics to attain a spatial resolution comparable to that of HST. 
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Fig. 1. — Radial velocity determination tests for 294 HST/STIS spectra of the calibration star, 
(a; left panel) Absolute difference It^diffl (in kms~^) between the measured radial velocity and 
true calibration star velocity as a function of the cross-correlation statistic Vcc- Open symbols 
indicate spectra with average signal-to-noise ratio S/N > 5.5 per pixel and closed symbols indicate 
spectra with S/N < 5.5. The vertical line indicates the cutoff = 2 that wc have used in the 
analysis, (b; right panel) Absolute difference |fdifr| (in kms^-"^) between the measured radial 
velocity and true calibration star velocity as a function of S/N. Open symbols indicate spectra 
with Tec > 2 and closed symbols indicate spectra with Vcc < 2. The vertical line indicates the cutoff 
S/N = 5.5 that we have used in the analysis. Both panels show a pronounced transition from 
unreliable measurements to reliable measurements, either with increasing or with increasing 
S/N. Spectra that have both Vcc > 2 and S/N > 5.5 (the open symbols to the right of the cut-off 
line in each panel) yield a reliable measurement of the stellar velocity. The horizontal line in both 
panels corresponds to 2.5 kms~^, which is the systematic uncertainty in the velocity calibration of 
the data (cf. Paper I). The axes in both panels are logarithmic. 
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Fig. 2. — Cross-correlation results for HST/STIS spectra of M15 stars. Spectra with an average 
signal-to-noisc ratio S/N > 5.5 per pixel were extracted for 131 stars. The ordinate of the figure 
shows the absolute difference It^diffl of the inferred stellar velocity and the systemic velocity of M15, 
divided by the central velocity dispersion of M15 (assumed to be 12kms-\ Gebhardt et al. 2000). 
This 'normalized velocity deviation' is shown as function of cross-correlation statistic r^c- A large 
value of the normalized velocity deviation generally indicates that the inferred stellar velocity is 
incorrect. There is a pronounced transition from unreliable measurements to reliable measurements 
with increasing Tqq. The cutoff at Vqq — 2 that we have used in our analysis (vertical line) excludes 
all measurements that are obviously unreliable. This confirms the analysis of the calibration star 
spectra presented in Figure 1. 
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B-V 

Fig. 3. — Color-magnitude diagram of the stars in M15, obtained from the photometric 
HST/WFPC2 imaging catalog presented in Paper I. The crosses indicate the location of the stars 
for which we have derived reliable radial velocities from the HST/STIS spectra. The two encircled 
crosses indicate potential RR Lyrae variables which were not included in the dynamical modeling. 
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Fig. 4. — Distribution of the number of M15 stars N with known Hne-of-sight velocities (with 
errors < 10 kms~^) as function of projected radius from the cluster center, (a; left panel) 
Binned histogram. Solid line: the HST/STIS velocity sample; Dashed line: the sample with 
ground-based measurements compiled by Gebhardt et al. (2000a). (b; right panel) Logarithm 
of the cumulative distribution (i.e., total number of stars within a given radius). Solid line: the 
HST/STIS velocity sample; Dashed line: the Gebhardt et al. (2000a) compilation; Heavy solid line: 
the combined HST/STIS and ground-based sample, corrected for double entries (i.e., stars present 
in both samples are counted only once). The HST/STIS data significantly increase the number of 
stars with known line-of-sight velocities close to the cluster center. 
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Fig. 5. — Positions of the stars with known hne-of-sight velocities in the central 4x4 arcsec of 
M15. Coordinates are (RA, DEC) with respect to the cluster center determined by Guhathakurta 
et al. (1996). Filled dots: stars with a velocity measurement from STIS; Circles: stars with a 
ground-based velocity measurement. The symbol size provides a measure of the brightness of the 
star; the circles were increased in size relative to the filled dots to avoid overlapping symbols. The 
STIS measurements cluster around a line with a position angle of 26.65°, which was the slit position 
angle used for the observations (see Paper I). 
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Fig. 6. — Comparison of the velocities derived with HST/STIS to those obtained from ground- 
based data, for those stars for which both exist. The excehent agreement shows that there are no 
unidentified systematic errors in the data. The two labeled stars are discussed in the text. 
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Fig. 7. — Velocity residuals, vhst — ground) for those stars that have both HST/STIS and ground- 
based measurements of their line-of-sight velocities, (a; top panel) Residuals as a function of 
stellar magnitude V. (b; bottom panel) Residuals as a function of broad-band color B — V . The 
weighted mean velocity residual is 0.2 it 1.4 kms~^. 
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Fig. 8. — Histogram of the differences Aujj between multiple velocity measurements of the same 
star. The histogram includes those stars for which multiple independent STIS spectra are available 
that individually satisfy the selection criteria with respect to blending, S/N and Vcc discussed in 
the text. The histogram is well fit by a Gaussian with a dispersion of 6 kms~^ (solid curve). The 
stellar velocities themselves can be estimated with an uncertainty that is factor 2 smaller than this 
(because the statistical uncertainty in the average of two measurements is a factor 2 smaller than 
the statistical uncertainty in their difference). 
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Fig. 9. — Radial profiles of projected kinematical quantities inferred from the combined HST/STIS 
and ground-based stellar line-of-sight velocity samples, (a; top left panel) Velocity dispersion 
cr; (b; top right panel) Rotation velocity V^-ot; (c; bottom left panel) Position angle PAkin of 
the kinematical major axis; (d; bottom right panel) The RMS projected linc-of-sight velocity 
fRMS = (^ot + where the angle brackets denote the average over a ring on the projected 

plane of the sky. Solid curves are best estimates; dotted curves define the 68.3% confidence bands. 
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Fig. 10. — Curves show the radial profile of ctrms with its associated uncertainty, as in Figure 9d. 

For comparison, data points show the absolute value of the difference Av = v^, — Vgys between the 
measured velocities of individual stars and the M15 systemic velocity. This allows some visual 
assessment of the radial trends in the data. The (Trms curves start at i? = 0.5", which is the 
average radius of the innermost 11 data points (11 is the adopted data window size in SROTATE). 
The radial distribution of the data points reflects not only the distribution of stars in M15, but 
also observational selection bias. 
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Fig. 11. — Data-model comparison for spherical dynamical models with an isotropic velocity 
distribution, a constant mass-to-light ratio T (in solar y-band units), and a central black hole 
of mass Mbh- (a; left panel) Likelihood contours as function of T and Mbh- The best-fit model 
has T = 1.6 ± 0.1 and Mbh = (3.2 ± 2.2) x 10^ Mq. Solid curves indicate the 1, 2, 3a and 4(7 
confidence regions in the two-dimensional (T, Mbh) plane. Dotted curves show likelihood contours 
inside the la region, (b; right panel) The RMS projected linc-of-sight velocity ctrms as a function 
of projected radius R. The heavy jagged curve surrounded by heavy dashed curves is the observed 
profile, as in Figure 9d. The smooth thin curves are the predictions for models with T = 1.6 and 
Mbh ranging from to 10 X W^Mq in steps of lO^M©. The dotted curve highlights the predictions 
for Mbh = 3 x 10^ M©, which is closest to the best-fit value Mbh = (3.2 ± 2.2) x 10^ M©. 
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Fig. 12. — Data-model comparison for spherical dynamical models with an isotropic velocity 
distribution, a mass-to-light ratio profile T(r) inferred from Fokker-Planck models (Dull et al. 1997), 
and a central black hole of mass Mbh- (a; left panel) The likelihood quantity A defined in 
equation (4) as function of Mbh- The minimum in A identifies the best fit black hole mass. 
Horizontal dashed lines indicate the 1 and 2(7 confidence regions, (b; right panel) The RMS 
projected line-of-sight velocity urms as a function of projected radius R. The heavy jagged curve 
surrounded by heavy dashed curves is the observed profile, as in Figure 9d. The smooth thin curves 
are the predictions for models with Mbh ranging from to 10 X 10^ Mq in steps of 10^ M©. The 
dotted curves highlight the predictions for Mbh = 4 x 10^ and Mbh = 5 x 10^ M©, which 
bracket the best fit value Mbh = (4.5 ± 2.1) x 10^ Mq. 
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Fig. 13. — Data-model comparison for spherical dynamical models with an anisotropic velocity 
distribution parameterized by equation (6), a constant mass-to-light ratio T, and no central black 
hole, (a; left panel) The likelihood quantity A defined in equation (4) as function of the central 
anisotropy Pq. At each /3o, the parameters a and T were varied to optimize the fit; (3oo was kept 
fixed at zero. The minimum in A identifies the best fit Pq. Horizontal dashed lines indicate the 1 
and 2a confidence regions, (b; right panel) The RMS projected line-of-sight velocity ctrms as 
a function of projected radius R. The heavy jagged curve surrounded by heavy dashed curves is 
the observed profile, as in Figure 9d. The smooth thin curves are the predictions for the best-fit 
models with /3o = 0, 0.2, 0.4, 0.65 and 0.8. The dotted curve is for /3o = 0.65, which provides the 
overall best fit. This indicates that considerable radial anisotropy must be invoked near the cluster 
center to explain the observations without invoking a central mass concentration. 
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Fig. 14.— The black hole mass -Mbh versus velocity dispersion <j . Solid points are measurements 
for various types of galaxies from the compilation of Tremaine et al. (2002). The solid line is the 
best linear fit, and the dashed lines show the Icr confidence band. The black hole mass estimate for 
M15 obtained from isotropic models, Mbh = (3.9 it 2.2) x 10^ (open point; this is the average 
of the masses inferred in Sections 5.3 and 5.4), fits perfectly onto the correlation. 
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Table 1. HST/STIS stellar velocity results 
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Note. — The ID number in column (1) corresponds to the number in the M15 stellar catalog presented in Paper I. The stellar positions [column (2) and 
(3)], V-band magnitudes [column (4)] and B — Y colors [column (5)] are taken from this catalog. The line-of-sight velocities inferred from the STIS spectra, 
and their uncertainties, arc listed in columns (6) and (7), respectively. The r^c cross-correlation statistic associated with the inferred velocity (see Section 2) 
is listed in column (8). The quantity / in column (9) measures the amount of blending in the spectrum (sec Section 3.3). The blending corrected velocities 
and their errors are shown in columns (10) and (11). Ground-based velocity determinations and their associated errors, from the compilation of Gebhardt 
at al. (2000a), are listed in columns (12) and (13). Ground-based velocity determinations with uncertainties in excess of lOkms""*^ were omitted (these were 
available for stars 5785 and 8760). 



